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ABSTRACT 

We present the Dynamic Eclipse Mapping method designed specifically to reconstruct 
the surface intensity patterns of non-radial stellar oscillations on components of eclips- 
ing binaries. The method needs a geometric model of the binary, accepts the light curve 
and the detected pulsation frequencies on input, and on output yields estimates of the 
pulsation patterns, in form of images - thus allowing a direct identification of the 
surface mode numbers (^, m). Since it has minimal modelling requirements and can 
operate on photometric observations in arbitrary wavelength bands. Dynamic Eclipse 
Mapping is well suited to analyze the wide-band time series collected by space obser- 
vatories. 

We have investigated the performance and the limitations of the method through 
extensive numerical tests on simulated data, in which almost all photometrically de- 
tectable modes with a latitudinal complexity £ — |m| ^ 4 were properly restored. The 
method is able by its nature to simultaneously reconstruct multimode pulsations from 
data covering a sufficient number of eclipses, as well as pulsations on components with 
tilted rotation axis of known direction. It can also be applied in principle to isolate 
the contribution of hidden modes from the light curve. 

Sensitivity tests show that moderate errors in the geometric parameters and the 
assumed limb darkening can be partially tolerated by the inversion, in the sense that 
the lower degree modes are still recoverable. Tidally induced or mutually resonant 
pulsations, however, are an obstacle that neither the eclipse mapping, nor any other 
inversion technique can ever surpass. 

We conclude that, with reasonable assumptions. Dynamic Eclipse Mapping could 
be a powerful tool for mode identification, especially in moderately close eclipsing 
binary systems, where the pulsating component is not seriously affected by tidal in- 
teractions so that the pulsations are intrinsic to them, and not a consequence of the 
binarity. 

Key words: stars:binaries:eclipsing - starsioscillations - asteroseismology - meth- 
ods:data analysis 



1 INTRODUCTION 

In recent years an increasing number of eclipsing binary sys- 
tems have been discovered to harbor pulsating components. 
The latest comprehensive catalog (Zhou 2010) contains 99 



such objects, discovered almost exclusively from the ground. 
The majority (70 in total) show delta Scuti type pulsations, 
followed by 8 sdB and 5 /3 Cephei type pulsators as the 



second and third most frequent types. Even more are ex- 
pected to be discovered from the current space-based mis- 
sions MOST, CoRoT, and Kepler. There are clear indica- 
tions that at least some of the modes are non-radial (as 
expected for such pulsators). For example, in systems like 
RZ Cas ( [Rodrfguez et al|2004[ ) or Y Cam ( |Kim et al|2002| ), 
the pulsations show amplitude and phase modulations dur- 
ing the eclipse phases, as a consequence of the symmetry 
violation in the surface flux integral due to the occupation. 
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This opens a new avenue of opportunities for asteroseis- 
mic investigations. As has been shown by numerous exam- 
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pies, there is a major difficulty in single stars that renders 
the mode identification a state of the art procedure: the 
inability to resolve their surface. The observables (fluxes, 
spectral line features) are weighted integrals of the local 
quantities over the visible stellar disc, and therefore show 
only weak dependence on the pulsation modes. All photo- 
metric and spectroscopic mode identification methods (for 
the most widely used ones, see Watson 1988j jBalona Ev-| 



2002 Briquet Aerts 2003; Zima 2006 ) must employ 



detailed models of the internal structure, atmosphere and 
pulsation, in order to overcome the problem of the low sen- 
sitivity. Accurate observational data and model parameters 
make a nearly unambiguous mode identification possible, 
with I and m determined to an accuracy of ±1 or better. 
Unfortunately, single stars rarely have well-determined pa- 
rameters. The more general approach of Doppler Imaging 
(Berdyugina et al. 2003 Kochukhov| |2004| , a remarkable 
technique for single stars, needs less sophisticated models 
because it aims at an image-like reconstruction of the sur- 
face patterns, and is therefore less sensitive to errors in the 
stellar parameters ( Kochukhov|2004 ^ ; however, it is only ap- 
plicable for rapid rotators, and the solutions still suffer from 
ambiguities (see, e.g., Berdyugina et al.]|2003 ). 

In contrast, a pulsating star in an eclipsing binary of- 
fers at least two advantages over the single star scenario. 
First, binarity enables a precise determination of the funda- 
mental stellar parameters. Second, the eclipses - the mutual 
occupations of the stars - implicitly provide a surface sam- 
pling: the shadow of one component literally sweeps across 
the surface of the other. This purely geometric phenomenon 
convolves the brightness distribution of the stellar surfaces 
into the variation of the integrated flux - that is, the light 
curve. Various inversion techniques may be used to recover 
the surface brightness structure from the light curve. They 
mainly differ in the amount of a priori assumptions about 
the surface pattern. The most common assumption, that 
the pulsations can be described by spherical harmonics, was 



employed by Gamarova et al. (2003), using the concept of 
spatial filtering (Nather & Robinson 1974p to identify non- 



radial modes in the interacting Algol- type binary RZ Gas. 
Although the modes could not be unambiguously identified 
- partly due to the complicating nature of the mass transfer 
between the components and partly because their approach 
was not suited very well to multiperiodic oscillations - their 
study was the ffist to demonstrate the potential of the ap- 
proach. Recently, Rodriguez et al (2010), based partially 
on the above method, made a preliminary mode identifica- 
tion for the 8 modes discovered in Y Gam, with a similar 
ambiguity in the mode numbers. 

More generally. Eclipse Mapping methods can be used 
to invert photometric time series into an instant image of the 
surface intensity distribution. When applied to high preci- 
sion photometric data with appropriate temporal resolution, 
they can discern far more detailed surface structure than the 
conventional methods. Eclipse Mapping techniques have al- 
ready been used to reconstruct static intensity patterns in a 
variety of eclipsing binaries (accretion discs in cataclysmic 
binaries, close binaries with spotted members, and even con- 
tact binaries) . An immediate opportunity is then to map the 
non-radial oscillation patterns in eclipsing binaries, which 
would allow a direct approach to the mode identification. As 
with Doppler Imaging, Eclipse Mapping needs only simple 



models with a few parameters, which, moreover, can be more 
easily determined in eclipsing binaries. Therefore, astero- 
seismology could in principle be made much easier for the 
pulsators in eclipsing binaries. In addition, the commonly 
employed approximation of non-radial pulsations with sin- 
gle spherical harmonics becomes questionable for moderate 
rotation speeds already (jAerts & Eyer 2000| , and is certainly 
invalid in rapid rotators (Lignieres et al. 2006"; 'Reese e t al.] 
[2006 J. Many single delt a Scuti stars are quite rapid rotators 
( Rodriguez et al.|2000 ), and rotation speeds of the same or- 
der are expected in binary systems with synchronized orbits. 
Obviously, Eclipse Mapping is a more realistic approach for 
such cases than fitting patterns of some analytic form. 

In the present paper we describe a variant of Eclipse 
Mapping method. Dynamic Eclipse Mapping^ designed to 
reconstruct the surface intensity patterns of nonradial pulsa- 
tions in eclipsing binary stars. We show that, under plausible 
circumstances, its application makes the mode identification 
possible in a large variety of eclipsing binary scenarios. 



2 DYNAMIC ECLIPSE MAPPING 

2.1 Eclipse Mapping of surface patterns 

The method of Eclipse Mapping (EM) was originally con- 
ceived to reconstruct the intensity distribution of radiation 
from accretion discs in eclipsing cataclysmic variables (GVs), 
with special emphasis on their radial temperature profile 



as a valuable diagnostic tool (Horne 1983 1985). Its over- 



whelming success in different GV scenarios (see Baptista 



2004| for a summary) has then inspired widespread applica- 
tions in other, less exotic eclipsing systems too, like mapping 
surface brightness inhomogeneities (spots) in late-type close 
binaries (e.g. XY U Ma,|Gollier Gameron|1997t or in W UMa 
stars (e.g. VW Gep, [Hendry et al.||1992| ). A simplified EM 
was also used to estimate the surface temperature map of 
an exoplanet occulted by its host star ( Knutson et al.|2006 ). 

The basic idea is that the eclipse acts as a surface sam- 
pler that convolves the surface brightness distribution into 
the light curve. Eclipse Mapping does the inverse, it per- 
forms a deconvolution of the light curve into a surface pat- 
tern. The latter is treated as an zma^^^omposed of pixels 
on an appropriate surface grid. 

No analytic assumptions are thus made on the surface 
brightness distribution, which allows the method to be max- 
imally free of stellar interior and atmospheric models, and 
makes it particularly suited to confront model predictions 
with observations in a maximally unbiased manner. All it 
requires is a proper modelling of the eclipses and some ba- 
sic atmosphere parameters. But because the convolution of 
a two-dimensional distribution into a one-dimensional time 
series implies a considerable loss of information, the inver- 
sion is ill-posed. The regularization is made by introducing 
additional, a priori information about the solution in the 
form of a regularization functional^ ^{f)^ which measures 
some desired property of the image /. The solution is ob- 
tained by maximizing this functional with respect to the 



^ Throughout the paper, the terms image, pattern and map are 
all used to refer to the same discretized surface brightness distri- 
bution. 
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elements of / and subject to the constraint of fitting the 
data at a prescribed level. The latter is usually measured by 
the chi-squared function. 

Classical Eclipse Mapping uses the information entropy 
of the image as the regularization function, 

<S(/, A) = - ^^^^ {h -Ak-fk Hh/Ak)) (1) 

(Home 1985| [Shore &; Joh nson 1980: Skilling"1988'), which 
measures the negative of the information content of the im- 
age vector f against that of a reference map A. The latter 
may be used to implement additional, user-defined prefer- 
ences, or may be just a uniform map, scaled so that its total 
flux equals that of /. The solution is then the image with 
the least structure (relative to the reference map) that can 
explain the observed data. 

Other regularization methods may be employed equally 
well; |Kaipio & SomersaloJ ( |2005 ) provide a thorough pre- 
sentation of possibilities, while Craig Brown| ( |1986| fo- 
cuses on inverse problems found specifically in astronomy. 
In particular, the choice of the regularization functional S 
depends on the type of the problem. For example, the in- 
formation entropy is a good choice for positive images with 
uncorrelated pixels. Physical maps are seldom uncorrelated; 
but the correlation can be easily accounted for via the refer- 
ence map (see Sec. 2.2 A\. Alternatively , the Tikhonov func- 
tional S{f) = II /V ('Tikhonov 1963b and its derivatives, 
which measure the smoothness of the image, are also used 



high surface degree ^, but, according to |Aerts, Christensen-| 
Dalsgaard & Kurtz (2010 Fig. 6.4.), photometrically de- 



in some cases, e.g., [Piskunov etH (1990). The optimization 
can be accomplished by standard methods. Most commonly 
it is transformed to a series of unconstrained optimizations 
by the method of Lagrange multipliers, where the multiplier 
plays the role of the regularization parameter, which is tuned 
between the optimizations until the desired level of data fit- 



ting is achieved. Skilling & Bryan ( 1984 ) give a sophisticated 
algorithm which accomplishes the two tasks in parallel, re- 
sulting in a robust and fast code. It is the algorithm that we 
have adopted in our implementation. 

2.2 Eclipse Mapping of pulsation patterns 

Previous applications of the Eclipse Mapping reconstructed 
surface structures that were static in time, or at least im- 
plicitly considered static during the data collection period. 
In principle only the underlying model needs to be changed 
in order to handle the time-dependent patterns. The time 
dependence, however, has implications on other aspects too. 



2.2.1 Assumptions 

The basic requirement of Eclipse Mapping is that the ge- 
ometric configuration of the binary must be known. Pul- 
sations, however, cause a periodic variation of the stellar 
shape, hence also in the local gravity and other atmospheric 
conditions. Their inclusion would require not only detailed 



atmospheric model (see Buta & Smith 1979 Townsend 



1997), but also the knowledge of the pulsation modes them- 



selves. Obviously, this is not feasible. Fortunately, for small 
amplitude pulsations these effects can be neglected. Indeed, 
for delta Scuti stars, the observed amplitudes in radial ve- 
locities, combined with typical pulsation periods, yield neg- 
ligible displacements compared to the stellar radius. Vari- 
ations of the surface normal could still be significant for 



tectable modes are limited by cancellation effect to ^ = 4. 
Although the eclipses break the symmetry of the disc inte- 
gration and may thus amplify some modes of even higher 
degree, those modes would still need to be detected outside 
the eclipses to have their frequency determined. Therefore, 
these variations can be safely neglected in our case, and the 
star is considered a rigid body of known shape. 

At the present stage we use spherical stars on circular 
orbits as the model of the binary. Limb darkening is taken 
into account (but its small variations with the local atmo- 
spheric conditions are again neglected) , but we make no fur- 
ther assumptions about the stellar atmosphere, using fluxes 
of arbitrary spectral range. 

For real, not-so detached systems, a proper account of 
the secondary's distorted shape may have to be made; al- 
though we note that its projection on the sky during the 
eclipses is still close to circular, so we limited ourselves to 
spherical secondary for the evaluation of the method's ca- 
pabilities. 

More stringent conditions are imposed on the pulsat- 
ing star, for which rotation is a complicating factor. Besides 
its physical influence on the pulsations, rotation has a sim- 
ple geometric effect that the frequency detected in the ob- 
server's frame will differ from the physical frequency in the 
co-rotating frame, according to the relation 



<^obs = <^surf + mQr> 



(2) 



where Q^rot is the angular rotation velocity, and m is the 
azimuthal order of the mode. Even if the rotation velocity 
were available, the physical frequency is unknown, because 
m is also unknown, being one of the parameters sought by 
our analysis. The time dependence of the patterns on the 
rotating stellar surface being thus unavailable, they cannot 
be mapped! What is known to us is the frequency of the 
patterns as seen on the visible, non-rotating stellar hemi- 
sphere. Therefore we have to map those patterns, and infer 
the whole-surface patterns from them. Obviously, this re- 
quires axial symmetry of both the stellar surface - in order 
for the shape of the disc to be constant in time - and of the 
pulsation patterns themselves - to assure that the surface 
and sky-projected patterns are equivalent, i.e., the pulsation 
amplitudes and phases seen on the visible stellar hemisphere 
correspond to those of the intrinsic pattern. 

Unfortunately, the above restriction excludes oscilla- 
tions of a tidally distorted star to be mapped with this tech- 
nique, because both the shape of the projected disc and the 
oscillation amplitudes do vary in time. Oblique pulsators, 
characteristic of roAp stars, would also be difficult to han- 
dle, the amplitudes being modulated by the stellar rotation 
(Kurtz 1982). In this respect, wider systems, where bina- 



rity and pulsation are a mere coincidence, are the preferred 
targets for the Dynamic Eclipse Mapping. 

2.2.2 Statement of the problem 

We may turn now to the mathematical formulation. For 
small amplitudes, each pulsation mode can be written as 
a sinusoidally oscillating perturbation to the equilibrium in- 
tensity distribution of the star, so the time-dependent (visi- 
ble) surface intensity pattern is the superposition of a static 



4 /. B. Biro and J. Nuspl 



equilibrium map, f'^^\ 
terns: 



and P sinusoidally oscillating pat- 



fit) = + Er=i [C*"' cosKt) + S*--) siii(w.t)] , (3) 
where uju is the frequency of mode z/, and the 'cosine' and 



sme maps 



A 



(^) 



cosF^^ and S*)^ 



(^) 



sinF, 



(^) 



have been used rather than the amplitude and initial phase 
maps (A, F) because they have the same units, of intensity. 

In the linear adiabatic approximation and for slow rota- 
tion, the pulsation patterns would be described in terms of 
spherical harmonics Yl^ . In the general approach, however, 
they are just pairs of images to be reconstructed. 

The dataset of integrated fluxes at M discrete moments 
is obtained by convolving the image /, composed of N pix- 
els, with an occultation kernel K, Suu M xN matrix, each row 
of which contains the pixels' contribution factors to the spe- 
cific datapoint. These factors are composed of the projected 
area of the visible pixel portion (including the foreshorten- 
ing factor), the limb darkening, and eventually other known 
factors incorporated into the model of the binary (proximity 
effects, for example). With the image ([3|, and denoting the 
times as t0((/) = 1 . . . M), the resulting synthetic dataset is 



(^) M^) 

\(t)k 



. (^) 



where we introduced the 'cosine' and 'sine' kernels 



(^) 

[c]cbk 



\ct)k 



= COS {(jJiytd)), 

= Kcf^k sin {(jj J, tcf,), 



(4) 



(5) 



obtainable from scaling the rows of the base kernel K with 
the corresponding cosine and sine time factors. With this no- 
tation. Equation Q can be written in a compact vectorial- 
tensorial form as: 



(6) 



Thus the known time dependence is transferred into the 
kernels, so that the model parameters are: i) one map for 
the static pattern, /^"^^ and ii) a pair of cosine and sine 
maps ( C^^\ S*^^-* ) for each pulsation mode, making together 
2P + 1 independent maps in total. 

The task is to estimate the 2P + 1 maps for a known 
kernel K, and a given set of frequencies cjz^ (z/ = 1 . . . P) as 
well as the light curve d. The kernel is provided by the model 
for the binary, while the frequencies and the light curve are 
the observed data. 



2.2.3 Regularization 

Because the pulsation patterns C and S are not strictly 
positive but may contain values of either sign, the entropy 
expression ([T]) is not valid for them (it can still be used for 
the static pattern f^'^^ though). For C and S we selected a 
simple quadratic regularization function, of form 



Sif,A) = -\\f-Af = -J2ih-Ak 



which also includes a reference map A. 



(7) 



This function is in fact a generalized Tikhonov func- 
tional, although we arrived at it through the statistical 



approach to inverse problems (Kaipio & Somersalo 2005), 



where the regularization function can be interpreted as the 
logarithm of an a priori probability distribution function 
(pdf) of the parameters. The above quadratic expression 
may be recognized as the logarithm of a joint Gaussian pdf, 
apart from a constant. The Gaussian is considered the most 
'generic' pdf (or most noncommittal with regard to missing 
information) when only the first two moments of the param- 
eter are available ( Jaynes|2003 Chap. 7 and 11). In our case, 
the first moment - the expectation value - is given by the 
reference map A, while the second moment - the 'spread' 
- plays the role of the regularization parameter. (Similarly, 
the entropy expression ([T]) is the logarithm of a Poissonian 
pdf, with mean pixel values given in the reference map A, 
e.g, i Skilling 1998 ) 

The entropy expression ([T]) for the equilibrium inten- 
sity map and the quadratic functionals ([7| for the pulsation 
patterns must then be simultaneously optimized, subject to 
a common constraint — Xaim- ^ P^i^ of maps belong- 
ing to one mode may be handled as a single entity, their 
regularization functions are simply summed up, so we have 
P + 1 objective functions to optimize, each on its subset 
of variables. This multiobjective optimization problem can 
be transformed with the simple weighting method to a sin- 
gle optimization of their weighted sum ^^^i "w^Su, with 
the weights determined by the desideratum that all maps 
should have about the same smoothness (as measured by 
the value of the regularization function). Because the objec- 
tives only interfere indirectly with each other via the com- 
mon data fitting constraint (they have disjunct subsets of 
variables) , the ideal point method can be employed for com- 
puting the weights i 



Liu et al. 2003| for the aforemen- 



tioned methods) . This involves performing P + 1 optimiza- 
tions first, with only one objective active at a time, but with 
all maps being fitted. A higher range of achieved values for 
an objective means more room to make the corresponding 
maps smoother, and translates to a higher weight in the final 
optimization run. 

We have found the Variable Chi Algorithm described 
by Baptista & Steiner ( 1993| quite useful in setting up a 



reliable criterion for data fitting. They introduced a so-called 
i?-statistics measuring the correlations of the neighboring 
residuals. It was shown that the R and statistics are 
proportional to each other, and setting a goal value for R 
as Raim 0.5 — 1 provides a more data-independent fitting 
criterion for - automatic noise scaling, in fact. 

The stopping criterion of the iterations was the same 
as in [Skilling Bryan^ (1984). After reaching a good fit- 
ting to the data, the iterations were continued until a TEST 
value, which measures the parallelism between the gradients 
of S{f) and C(/, d), decreased below a certain value in all 
subsets. Being half the sine of the angle between the two in- 
volved vectors (in terms of the scalar product), TEST takes 
values from to 0.5. The algorithm could routinely reach 
below TEST=0.01, which we chose as the stopping value. 



2.2.4 The reference map 

The role of the reference map A appearing in the regular- 
ization functionals is to allow the introduction of additional 



Photometric mode identification 1. Dynamic Eclipse Mapping 5 



user preferences about the solution. In highly ill-posed in- 
verse problems like Eclipse Mapping, where every bit of a 
priori knowledge is important, it may have a determining 
role in obtaining a proper solution. 

The regularization functionals have their global maxi- 
mum at the location of the reference map, which therefore 
is the default solution in the absence of observational con- 
straints. Regular solutions will also be as close to it as al- 
lowed by the constraints. Updating the reference map from 
the instantaneous solution during the iterative solving pro- 
cedure is a common technique to make it control the im- 
age property to be measured by the regularization function. 



which is hence optimized in the solution (see Horne||l985| 
and |Bobinger et aL]|1999| for examples in accretion disc re- 
constructions) . Alternatively, it may be set to a uniform map 
if no additional preferences exist; in this case we obtain a 
'most uniform', but spatially uncorrelated, solution: permut- 
ing the pixels will not change the value of the regularization 
function at all. This is better than nothing; but most prob- 
lems do have some symmetry, and employing it improves the 
solution. 

The pulsation patterns of a rotationally symmetric star, 
when described in a spherical coordinate system tied to the 
rotation axis, obey a kind of axial symmetry, in that the 
local amplitude only depends on the latitude, while the lo- 
cal phase only varies with the longitude. This holds not 
only for spherical harmonics, but for all pulsations of tidally 
izndistorted stars - including fast rotators, for example. 
Consequently all the 'cosine' and 'sine' maps are expected 
to have the form C{0,ip) = A{0) cos F {(p) , and S{0,(p) = 
A(0) sin F ((f) , where and cp are the co-latitude and lon- 
gitude, respectively. On a uniform spherical grid, therefore, 
we have a discrete number of amplitudes Ak = A{Ok) {k = 
1 . . . Ne) and initial phases Fi — F((pi) (/ = !... N^), that 
completely describe the maps C and S on that same grid: 



Cki = Ak cos Fi, 
Ski = Ak sin Fi. 



(8) 



According to this a priori expectation, we set up the 
following updating scheme for the reference maps. The al- 
gorithm starts with computing flat-valued images by linear 
least-squares fitting to the data. The reference maps are ini- 
tialized with the same images. After the end of each iteration 
the latest (C, S) solutions are interpolated onto a uniform 
spherical grid around the rotation axis, of appropriate reso- 
lution Ne X N^. The non- linear model ^ with the param- 
eters Ai . . . Anq, Fi . . . Fn^ is fitted to the maps, by a stan- 
dard Levenberg-Marquardt algorithm. Spatial correlations 
are then taken into account by smoothing the resulting dis- 
cretized amplitude and phase profiles with a Gaussian of a 
user-supplied angular correlation length. From the smoothed 
profiles, new (C, S) maps are computed and interpolated 
back to the reconstruction grid, giving the reference maps 
for the next iteration. For the equilibrium map f^'^^ only 
a surface Gaussian smearing is applied; this takes spatial 
correlations into account (i.e., to yield a smooth image). 

A few safety measures were needed for proper operation 
of the above procedure. These include a careful estimation 
of starting parameter values for the non-linear model fitting 
algorithm, as well as making the smoothing wraparound for 
the phase profiles (running along the longitude) and mir- 



rored at the poles for the amplitude profiles, in order to 
ensure continuity over the stellar surface. 

The same procedure is used when amplitude and phase 
profiles are derived from the final solution for the purpose 
of mode identification, with the difference that the profiles 
are not smoothed. 

The above scheme only works for a known rotation axis. 
Stars in eclipsing binaries are generally assumed to have the 
rotation axis perpendicular to the orbital plane. This as- 
sumption is expected to hold for close binaries, where tidal 
interactions tend to bring the spin axes parallel to the or- 
bital axis, in addition to circularizing the orbits. For wider 
systems, which are of primary interest for us (Sec. 2.2.1), it 



may not be true, as demonstrated by the recent discovery 
of very oblique rotation axes (being almost in the orbital 
plane) in the eccentric system DI Her ( Albrecht et al.|2009 ). 
If the rotation axis is unknown, then only the assumption of 
smoothness can be used, leading to a 'smoothest' solution, 
which will for certain be inferior to the regular solution, 
and may be insufficient for mode identification. Assuming a 
wrong axis would do even worse, though. 

We note that the seemingly ad-hoc manipulation of 
the reference map may be justified within the frame of the 
Bayesian approach: each iteration is a prior-to-posterior pro- 
cessing step, and the reference map updating only prepares 
the prior for the next iteration from the posterior of the 
last iteration, by propagating only those structures that are 
consistent with the a priori expectations, and ignoring ev- 
erything else. 



2.3 The role of the eclipse geometry 

The amount of recoverable surface structure information 
depends primarily on how the surface is sampled by the 
eclipses. Uneclipsed regions, for example, cannot be recon- 
structed; but we cannot ignore them either because they 
contribute to the integrated fiux. Therefore pixels within 
these regions are fitted, but not 'regularized'. Being spa- 
tially unresolvable, the uneclipsed region is replaced by a 
single parameter, a virtual pixel, with its contribution factor 
equivalent to that of the whole region. The value of this pixel 
is simply copied to the virtual pixel of the corresponding ref- 
erence map during the updating scheme, therefore it gives 
zero net contribution to the regularization function, but still 
counts to the goodness-of-fit. This technique also prevents 
an unnecessary proliferation of the parameters. Each map 
has its own virtual pixel of this kind; the one pertaining to 
the static intensity map takes in addition the duty of fitting 
any other uneclipsed fiux ('third light') in the system. 

Although systems with total eclipses seem at first sight 
the most favorable scenery for EM because they sample the 
whole stellar surface, the actual situation is more subtle. On 
one hand, there is no useful data from the totality phases 
of the light curve, so the amount of information may be less 
than with partial eclipses. On the other hand, higher inclina- 
tions - for which the eclipses become total - do not necessar- 
ily imply a better reconstruction. As the system approaches 
the edge-on configuration, the eclipses start to sample mul- 
tiple areas in exactly the same way. This latter is illustrated 
in Fig. [l] The ingress and egress arcs, corresponding to the 
companion's projected limb on the stellar disc at various or- 
bital phases, form a 'sampling grid'. Each pixel of the grid 
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of such an optimal range requires that 



which 



Figure 1. The sampling grids of the primary's sky-projected 
disc for a binary system with fractional radii Ri = 0.195 and 
R2 = 0.265, seen under inclinations of 80° (left panel) and 85° 
(right panel), respectively. The secondary sweeps across the pri- 
mary in the horizontal direction. The white upper areas are never 
eclipsed. Each shade of the grey-coloured pixels marks one spe- 
cific chain of equivalent pixel pairs on the right panel. One such 
pair is highlighted with thick black lines. Both panels also show 
a pixel with no equivalent pair. 



is eclipsed at one phase interval and reappears at another, 
contributing to exactly two data points of the differential 
light curve, interpreted in terms of flux changes from phase 
to phas^ The left panel shows a configuration in which all 
the ingress and egress arcs intersect each other in at most 
one point. The right panel in turn presents a configuration 
in which most arcs have two intersection points, so there are 
pairs of pixels that are eclipsed simultaneously and reappear 
simultaneously. No inversion method will be able to separate 
their contribution from each other. 

It is easy to see that these equivalent pixel pairs are ar- 
ranged symmetrically with respect to the trajectory of the 
secondary's centre, projected on the sky. If the trajectory 
intersects the primary's disc, then there will be a horizontal 
band similar to that of the uneclipsed region, full of equiv- 
alent pixel pairs, and therefore it is an ambiguous region, 
for which the reconstruction is likely to be distorted. The 
remaining area in between represents the trusted region, for 
which the unique sampling makes a trustworthy reconstruc- 
tion possible. 

The vertical extents of the limiting regions can be eas- 
ily computed for spherical stars. They are also approxi- 
mately valid for real stellar shapes, if polar radii are used. 
Thus, for stellar radii Ri and R2, binary separation a, 
and inclination i, the uneclipsed region has a fractional 
height hv = 1/2 — (R2 — acosi) / (2Ri), while for the am- 
biguous region it is /ia = 1 — acosi/Ri, all in units of 
the stellar diameter. For inclinations larger than imin — 
arccos((i?2 — Ri)/a), the eclipses are total (formally hv<0). 
The ambiguous region disappears for inclinations smaller 
than imax = arccos(i?i/a). Inclinations between these two 
limits correspond to the best cases, in which the whole stel- 
lar disc is unambiguously sampled. However, the existence 



^ In fact, reconstructing on the sampling grid would be most 
welcome, because its elements hold the discrete 'packets' that are 
coded into the light curve. Unfortunately, the pixel areas vary 
significantly on such a grid; moreover, the grid depends on the 
sampling of the observed dataset, which may change from eclipse 
to eclipse - although the latter would be easily overcome by re- 
sampling the data on a uniform temporal grid. 



implies R2 ^ 2Ri; that is, the secondary must be at least 
twice as large as the primary. Although such systems may 
exist, in most cases the secondary is not so much larger. In 
addition, not all systems show total eclipses. The reconstruc- 
tion therefore will always be compromised to some extent by 
the limits of the eclipse geometry. But since partial eclipses 
may not be as bad compared to total eclipses as it may seem 
at first, any inclination for which the whole eclipsed region 
is uniquely sampled is close to optimal. 

If the eclipses are central (i = 90°), Eclipse Mapping 
will only reconstruct modes that are symmetric with re- 
spect to the orbital plane. In an aligned rotator, asymmet- 
ric modes (I — m — odd) are subject to complete cancel- 



lation (Chadid et al. 2001), which also persists during the 
eclipses due to their symmetry, so only symmetric modes 
are detected anyway. However, in an oblique rotator (not 
pulsator!), non-symmetric modes will be sampled in a sym- 
metric way, and therefore will not be properly reconstructed. 
Such edge-on systems are expected to be rare, though. 

From Fig. ^ it is also clear that the resolution of the 
sampling decreases in the vertical direction, so the southern 
regions (from the observer's point of view) will be recon- 
structed with poorer quality. There is no such intrinsically 
uneven sampling in the horizontal direction. Hence, for an 
aligned rotator, the phase profile will be reconstructed rea- 
sonably well, which allows a reliable identification of m; but 
the amplitude profile will suffer from imperfections in the 
southern parts. Fortunately, for the rotationally symmetric 
pulsations appropriate for the method it suffices to know 
the pattern on one half of the stellar hemisphere in order to 
identify the modes, so the eclipse geometry might not have 
a dramatic impact on the mode identification. On the other 
hand, an oblique rotator has more subtle implications, and 
the above arguments are not applicable; both profiles will 
be affected to some extent by the uneven vertical sampling. 

To investigate how the limitations discussed above af- 
fect the reconstructions, we have selected four system con- 
figurations for testing (two distinct systems viewed under 
two different inclination angles each), listed in Table [l] and 
representing various levels of optimality, as indicated by the 
heights of the key regions. System 1 is an ideal system, with 
the secondary more than twice as large as the primary and 
an inclination which enables unambiguous sampling of the 
whole surface and makes the eclipses total. System 2 is the 
same system seen at a higher inclination, with an ambigu- 
ous region. System 3 features both an uneclipsed and an 
ambiguous region, while System 4 has almost total eclipses 
but an ambiguous region of considerable area. 



3 TESTING THE DYNAMIC ECLIPSE 
MAPPING 

In order to asses the performance of the method, we have 
subjected it to extensive numerical testing on artificially 
generated data. The model uses a simple binary geome- 
try, with rigid spherical stars on a circular orbit. Likewise, 
we employ a simple model for the stellar atmosphere: bolo- 
metric limb darkening, with coefficients taken from | Claret | 
( |2001 ) for an average delta Scuti pulsator with T — 7500 K, 
\ogg = 2.5, solar composition, and no turbulent velocity. In 
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Table 1. Binary system configurations selected for reconstruction 
tests. 



System 


Ri 


R2 


i 




Ha 


htr. 


1 


0.153 


0.352 


79.4 








1 


2 


0.153 


0.352 


83.6 





0.24 


0.76 


3 


0.195 


0.265 


80.7 


0.23 


0.17 


0.60 


4 


0.195 


0.265 


85.1 


0.04 


0.56 


0.40 



R2 are stellar radii in units of the binary separation, hy and 
/lA cire the fractional heights of the uneclipsed and ambiguous 
regions respectively; the remaining amount is the trusted region 
with height htr. • 



most runs, however, only the linear limb darkening law was 
used, with the coefficient xi — 0.56 drawn from the linear 
part of the polynomial relation. 

For strictly spherical stars the mutual occultations can 
be analytically computed. For maximum precision, each 
pixel is composed of two flat, triangular tiles, where all the 
quantities (including partial visibilities due to the occulta- 
tions) are computed separately and then summarized for the 
whole pixel. 

To avoid the so-called inverse crimes, a situation which 
occurs when the data generation and the inversion are made 
with exactly the same setup, thus ignoring the modeling 
errors and leading to an over-optimistic assessment of the 
method's performance (Kaipio Somersalo| |20()5 ), we used 
two different grids for the direct and inverse parts. For gen- 
erating the artificial light curves, a uniform spherical grid 
aligned with the rotation axis was set up on the stellar sur- 
face, with an exaggerated resolution of 240 x 120 grid ele- 
ments. In turn, the reconstruction was done on an adaptive 
polar grid, applied on the visible stellar hemisphere, with its 
z axis pointing towards the observer (as discussed in Sec- 
tion |2.2.1| the reconstructions must be done on the fixed 
visible hemisphere). The grid was chosen so that its projec- 



tion on the sky consists of concentric rings of equal widths, 
each ring being divided into a number of identic pixels de- 
termined by the condition that their contribution factors to 
the integrated flux, wis., pixel area x cos 7 x L(cos7) (where 
7 is the aspect angle and L is the limb darkening) , should be 
as uniform as possible across the grid. With this particular 
choice, the average linear dimension of the pixels is also of 
the same order, and can be easily matched to the average 



resolution of the sampling grid (Sec. 2.3), ensuring that the 
quality of the reconstruction is not limited by the pixeliza- 
tion. In addition, the approximate constancy of the above 
contribution factors means that the reconstruction of surface 
intensities and projected fluxes of the pixels are practically 
equivalent, and in fact the latter case, being technically sim- 
pler, has been implemented in our algorithm. 

Throughout the reconstructions we used a grid com- 
posed of 30 such rings, summing up to ^ 3000 pixels cover- 
ing the stellar disc. 

As demonstrated earlier inlBiro & Nuspl (2005), in prin- 



ciple it is possible to reconstruct the static equilibrium in- 
tensity map along with the pulsations. In practice, however, 
the appropriate manner is to account for any static flux com- 
ponent during the binary model fitting; otherwise it would 
probably spoil the binary parameters. (The only purpose 
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Figure 2. Illustration of the mode identification procedure. Pan- 
els (a) and (b): the 'cosine' and 'sine' maps of the reconstructed 
projected stellar disc for a single mode (here £ = 3 and m = 1); 
crosses mark the position of the rotation axis. Panels (c) and (d): 
the equivalent Mercator maps in the spherical system of the ro- 
tation axis. Isocontour lines are used to enhance the structure 
on the grayscale maps; thin black and white isocontours mark 
positive and negative values, respectively. The zero isocontours, 
corresponding to the node lines, are drawn with thick black lines. 
The white "polar caps" are due to the uneclipsed region. Pan- 
els (e) and (f): the amplitude and phase diagrams fitted to the 
maps (c) and (d). Thick solid lines are the profiles of the input 
model, those of the reconstruction appear as thin solid lines, with 
the error estimates of the individual points drawn with vertical 
lines. Dashed vertical lines in panel (f) show the longitude range 
corresponding approximately to the visible stellar hemisphere. 



of re-mapping it with EM would be during an iterative re- 
finement of the parameters.) Therefore we consider it as a 
nuisance factor and leave it out from the current analysis. 

We generated evenly sampled synthetic light curves cov- 
ering a range slightly larger than the eclipses. With typical 
pulsation frequencies being 10 — 100 times the orbital fre- 
quency, we chose a sampling interval of 0.00163, in units of 
orbital phase. One cycle of the fastest pulsation is thus cov- 
ered by about 6 data points. With typical system configura- 
tions this amounts to about 150 data points per eclipse. Up 
to 20 eclipses from successive orbital cycles were involved. In 
most cases, though, a much smaller number (sometimes even 
a single eclipse) already yielded a successful reconstruction. 

The observational errors were modelled by adding arti- 
ficial Gaussian noise of a specified level to the synthesized 
light curve. In this study we use a detection signal to noise 
ratio (snr), the signal being not the total flux, but rather the 
semi-amplitude of the weakest pulsation mode, as measured 
in the total flux outside of the eclipses, where it is free of dis- 
tortions. The noise is 3cr, as usual (a is the parameter of the 
Gaussian noise). Typically we used values of snr = 2 — 10. 

3.1 Mode identification 

For each pulsation mode. Dynamic Eclipse Mapping recon- 
structs a pair of intensity patterns on the adaptive grid ap- 
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plied on the visible stellar hemisphere. They are shown in 
panels (a) and (b) of Fig. |2] with the addition of limb dark- 
ening for visualization purposes - to look exactly as a close 
observer would see them in the absence of the other modes, 
at a given moment t = and a quarter of period later, 
respectively. These data are meaningful after transformed 
into the spherical coordinate system of the stellar rotation 
axis (be it given or assumed), where they can be visualized 
in the usual form of Mercator maps, as illustrated in pan- 
els (c) and (d). These maps have the special property that 
their half corresponding to the invisible stellar hemisphere 
is empty Recall that the reconstructions are made on a 
steady hemisphere, the rotation being redeemed by symme- 
try assumptions. Eventual uneclipsed regions, lacking any 
spatial information on the pulsation pattern, will appear 
also as empty "polar caps" near the northern pole of the or- 
bital axis, their fluxes being accounted for by virtual pixels. 

In most cases the pulsation mode can already be guessed 
by visual inspection of the node lines (isocontours of zero 
value). However, their identification can be put on more 
quantitative grounds, by fitting the reference map model ^ 
to the final solution in order to get amplitude and initial- 
phase profiles, shown in panels (e) and (f). (Due to the miss- 
ing base intensity map, the pulsation patterns have no abso- 
lute scale, and that is why only the zero level is marked on 
the amplitude profile diagrams.) The slope of the phase pro- 
file in the central parts of the stellar disc (|(/?| ^ 60°, say), 
rounded to the nearest integer, gives the azimuthal mode 
number m, and at the same time imposes a lower limit on I. 
Then checking the number and positions of the nodal points 
(roots) of the amplitude profile within the trusted region, 
and comparing them with those of the possible associated 
Legendre polynomials PP{cosO), gives the most probable 
i ^ \m\ value. We followed this simple procedure in deter- 
mining the mode numbers in each case. 

The underlying non-linear fitting algorithm used for fit- 



ting the profiles (Levenberg-Marquardt, Sec. 2.2.4) also fur- 
nishes error estimates for the amplitude and phase profiles. 
However, we found that sample statistics on columns and 
rows of the polar grid used for their computation give more 
reliable estimates, therefore we used the latter for the errors. 



3.2 Single modes 

We performed reconstructions of single pulsation modes up 
to ^ = 4 on an aligned rotator in all the four systems listed 
in Table [l] A frequency of 63.1234 cycles per orbital pe- 
riod and randomly chosen initial phases were used. Artificial 
data covered the phases between ±0.11 around 1 to 5 suc- 
cessive eclipses, and had a detection signal to noise ratio of 
10. Noiseless sample modulation light curves in System 1 are 
shown in Fig. [3] The orbital phases affected by the eclipse 
are approximately between ±0.082; as it can be seen, the 
eclipses are total for |(/)orb| ^ 0.02. 

In most cases, data from a single eclipse were enough 
for a reliable reconstruction, although a few of the modes 
- notably the zonal modes with i 2 - required at least 5 



^ Although the invisible hemisphere could be interpolated from 
the visible hemisphere once the mode was identified, that would 
be confusing and unrealistic. 
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Figure 3. The modulation of the pulsations by the eclipses in 
System 1: single-eclipse excerpts from noiseless synthetic data. 
Panels (a) and (b) summarize the zonal and sectoral modes; the 
tesseral modes are shown in panels (c) and (d). The mode num- 
bers are shown on the left of the light curves. The fittings achieved 
by Eclipse Mapping are drawn with solid lines. 



Fig. [4] shows the reconstructions for a selection of 



eclipses. Therefore all the cases were investigated with this 
eclipse coverage 

modes. □ 

For completeness, we repeated some of the reconstruc- 
tions using other frequencies, about three times lower and 
higher than our adopted value, respectively. The results were 
the same in all cases, meaning that the pulsation frequency 
does not play a significant role in the reconstruction, as long 
as the pulsation is adequately sampled. 

In general, the features of the reconstructions are in 
agreement with the limitations discussed in Sec. |2.3| The 
phase profiles are uniformly restored over the central part 
of the stellar disc in almost all cases, as clearly demonstrated 
by the amplitude and phase diagrams of Fig. |4] for longi- 
tudes within ±60°, almost all the phase profiles run together 
with the model. There is only one exception: mode (4, 1) in 
Systems 2 and 3. For the amplitude profiles, the reliability of 
the reconstruction varies with the latitude as predicted; the 
northern half is far better reproduced than the southern. In 
addition, a peculiar reversal of the patterns may be observed 
near the north pole, leading to an 'overshooting' in the am- 
plitude profiles, which causes false nodes near 1^ = 0°. This 
phenomenon appears when the input models have a node at 
the pole; it is absent from the zonal modes. Fortunately the 
photometrically detectable modes with ^ ^ 4 do not have 
nodal lines so close to the polar regions, so this artefact 
should cause no confusion. 

For each investigated system there were a few cases 
for which the reconstruction yielded false mode numbers, 
at least with the default setup of 5 eclipses. We list them in 
Table[2] While System 1 is indeed the best of all as expected, 
the differences are not large. It seems that neither the large 
uneclipsed regions (as in System 3) nor the ambiguous re- 

^ The full set of reconstructions, including modes with negative 
m, is available in the online material. 
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Figure 4. A selection of reconstructed single modes in Systems 1-4, aligned rotator case. Each column contains data for a particular 
mode, marked at the top of the columns, and ordered from left to right in increasing complexity measure of their vertical structure, 
i — \m\ (see text). The 'cosine' maps are shown for each case, with contour lines overlaid, as in Fig.|2] Row (a) shows the input model, 
rows (b) to (e) contain the reconstructions in Systems 1 to 4. Cases with unsuccessful mode recovery are labelled by the misidentified 
mode numbers. The lowest two rows (f) and (g) show comprehensive amplitude and phase profiles for all systems. In these diagrams, 
thick dashed lines correspond to the model profiles, while the reconstructed ones are drawn with thin solid, long dashed, short dashed 
and dot-dashed lines for Systems 1-4, respectively. 



gion (as in System 4) cause any dramatic degradation of the 
profiles. Our conjecture is that the specific reference map 
updating scheme has a significant role in this rather nice 
behavior. 

Looking at Table [2] one may observe a decreasing qual- 
ity of the reconstruction with increasing i — \m\. This dif- 
ference is the number of latitudinal nodal lines in the ampli- 
tude profiles (not counting the polar nodes), and is therefore 
a rough measure of the pattern complexity in the latitudi- 
nal direction. However, even the most complex mode, (4, 0), 
could be recovered in 3 out of 4 systems. 

In general, it can be concluded that the success of mode 
identification depends primarily on the mode's complexity 
in the direction perpendicular to the sampling direction of 



the eclipses. That said, in the case of a primary with the 
rotation axis tilted by 90° sideways (that is, in the plane 
of the sky), for instance, zonal modes would be the most 
favorable cases, while the algorithm would have difficulty in 
restoring mode (4,4), for example. 

3.3 Towards reality: multimode pulsations and 
oblique rotators 

Single mode nonradial pulsation is of course an ideal case. 
Real pulsators have at least a dozen simultaneous modes. 
Given the frequencies of the detected modes, however, 
Eclipse Mapping is able in principle to separate each mode's 
contribution from the others. Fig. [5] presents a test case 
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Table 2. List of misidentified modes in the test systems 





dv^ cri "n 1 
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(4,0) 


(3,0) 
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(3,0) 


(4,0) 




(3,1) 


(2,1) 




(4,1) 


(4,2) or (4,3) 




(4,2) 


(3,2) 
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(1,0) 


(4,0) 




(4,1) 


(4,3) 
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(3,0) 


(4,0) 




(4,0) 


(3, 0) or (4, 0) 



Table 3. Parameters of the simulated three-mode pulsational 
case. The amphtudes are in arbitrary units and refer to the semi- 
amphtude in the integrated flux outside the echpses. 



(1,0) (3,1) (2,2) 




Figure 6. Echpse maps of a triple-mode pulsation of an oblique 
rotator in systems 1 and 3 (rows 'a' and 'b', respectively). Mode 
numbers are indicated at the top of each column. Shown are the 
'cosine' maps of the projected stellar disc. The symbols mark 
the approximate position of the rotation axis. 



Mode Frequency Amplitude Initial phase 

[^orb] [°] 



(1.0) 59.153517 1.387 143.56 

(3.1) 61.547029 1.000 65.44 

(2.2) 65.787702 1.839 12.05 



Table 4. Frequency sets for multimode reconstructions. 



Set uji UJ2 CJ3 



1 59.153517 61.547029 65.787702 

2 59.0 61.0 63.0 

3 59.153517 61.153517 63.153517 



of three simultaneous modes with parameters listed in Ta- 
ble |3] and reconstructed on an aligned rotator in Systems 1 
and 3. (Sky-projected images are shown in order to allow a 
comparison with the oblique rotator case below). The fre- 
quencies were intentionally chosen close to each other, to see 
how well their modes can be separated by the method. We 
also decreased the signal to noise ratio to 5. Accordingly, 10 
eclipse cycles had to be included for a good reconstruction. 

Similarly, in real systems the axis of rotation is not nec- 
essarily aligned with the orbital axis, especially for wider 




Figure 5. Eclipse maps of a triple-mode pulsation of an aligned 
rotator in systems 1 and 3 (rows 'a' and 'b', respectively). Mode 
numbers are indicated at the top of each column. Shown are the 
'cosine' maps of the projected stellar disc. The x symbols mark 
the approximate position of the rotation axis. 



binaries with weak tidal interactions. In this case, how- 
ever, the direction of the rotation axis must be known a 
priori. Assuming that this is the case, the specific refer- 
ence map scheme provides a reliable reconstruction for tilted 
axis as well. Fig. [6] shows the repetition of the triple-mode 
pulsation case for a rotation axis tilted with Eulerian an- 
gles (0, ^,'0) = (41°,63°,0°) - that is, rotate the meridian 
around the initial axis by 41 degrees, then tilt the axis in 
the new meridional plane 63 degrees towards the equator. 
The third rotation around the new axis would be equivalent 
to a shift in the initial phases of the pulsation patterns, and 
was therefore omitted. As a side effect, the longitudes cor- 
responding to the visible disc shift from [—90, 90] to about 
[-180,0]. 

Based on the topology of the node lines, it can be seen 
that the mode numbers (1,0), (3,1) and (2,2) are success- 
fully recovered in both cases. This is also confirmed by the 
amplitude and phase profiles (available in the online mate- 
rial) . 



3.4 The role of multiple eclipses 

If the frequencies are not in resonance either with the orbital 
motion (as would be the case for tidally induced oscillations) 
or with each other (as would occur for mode interaction), 
then, for every mode, the surface patterns at a given or- 
bital phase will differ from eclipse to eclipse; therefore every 
new observed eclipse will present independent pieces of in- 
formation to the algorithm, leading to a better separation of 
the modes, thus improving the quality of the reconstruction. 
The increasing amount of independent information will also 
compensate for the lower quality of a noisy dataset. 
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Conversely, for tidally excited pulsations, where the fre- 
quencies are multiples of the orbital frequency, the flux mod- 
ulation repeats itself exactly from orbit to orbit, making 
the additional eclipses largely redundant. A similar scenario 
is the frequency splitting of non-radial modes of a syn- 
chronously rotating primary, where the differences of the 
frequencies are multiples of the orbital frequency. Although 
the rotational splitting is not strictly uniform due to addi- 



Table 5. List of three equivalent binary parameters giving similar 
primary eclipse profiles. 



tional second- and higher order factors (Soufi et al. 1998 
Goupil|[2000 ), nevertheless they will be quite close to reso- 



nance, and pose the same problem, this time the repetition 
of the relative phase differences between the split modes. In 
these cases, multiple eclipses would not be able to improve 
the reconstruction above a certain limit. Their only benefit 
will be the improvement of the observational signal-to-noise 
ratio. 

To investigate the above possibilities, we repeated the 
oblique rotator case for System 1 (to keep the complexity 
level reached so far) with two other frequency sets, simu- 
lating the two resonance possibilities discussed above. They 
are listed in Table [4] with the original one in the first line. 

The results are compared in terms of the reconstructed 
profiles in Fig.|7| In the non-resonant case, two eclipses still 
give poor results, in particular for the first mode (row (a)), 
but raising the number of eclipses to 10 (row (b)) provides 
enough information for a proper mode identification. In con- 
trast, the results for the resonant cases, shown in rows (c) 
and (d) for the same number of 10 eclipses, demonstrate 
that the improvement of the reconstruction expected from 
the inclusion of additional eclipses is prevented by the reso- 
nances. In our particular case, the fitted profiles lead to the 
estimates (3, 1), (3, 1), and (4, 2), instead of the true figures 
(1,0), (3,1), and (2,2). 

Tidally induced oscillations and synchronous rotation 
occur primarily in close binaries. The mapping in such sys- 
tems is already complicated by other proximity factors (tidal 
distortion effects, mass transfer), and the limitation dis- 
cussed above is one more reason why wider binary systems 
are more preferred by Dynamic Eclipse Mapping. It should 
be emphasized, though, that once the effects of interaction 
are properly accounted for, then the close binaries with non- 
resonant pulsations become more attractive, because the 
components are more likely to contain aligned rotators - 
that is, with a known direction of the spin axis. 

3.5 Errors in the modelling parameters 

All the tests made so far assumed the true values for the 
physical parameters. Since EM by its nature is sensitive to 
the accuracy of the data, it is sensible to ask what impact 
do inexact model parameters have on the results. Parame- 
ters that may affect the reconstructions are: the frequencies 
of the pulsation, the geometric parameters (stellar radii, in- 
clination, and the direction of the spin axis), and the limb 
darkening. 

Recent space-based, continuous, long-term observations 
allow a very precise determination of the pulsation frequen- 
cies, therefore serious errors in the frequencies shouldn't 
happen these days. 

Likewise, the geometric parameters can in general be 
determined to a high accuracy in eclipsing binaries. How- 
ever, there are well-known correlations between the stellar 
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[sep] 


[sep] 


[deg] 
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0.153 


0.352 


79.4 


2 


0.160 


0.325 


82.4 


3 


0.173 


0.302 


85.4 



radii and the inclination, in that about the same eclipse pro- 
file may be reproduced by different combinations of their 
values. With poor or missing spectroscopic data, more pa- 
rameters have to be 'guessed', and the correlations might 
lead to biased parameters. The problem of incorrect or un- 
known direction of the rotation axis is more complex, and 
is left for a subsequent paper. 

To investigate the effects of correlated deviations from 
the true values, we created three sets of binary parameters 
by shifting the inclination by 3° in both directions from a 
central value, and modifying the stellar radii accordingly by 
a trial and error method, to get about the same primary 
eclipse profile. The parameters are listed in Table [5] The 
limb darkening was kept fixed. We then generated data for 
a triple mode pulsation in the middle system with i — 82.4°. 
We used the same pulsation model as in Table js] with the 
difference that the third mode (2, 2) was replaced by (4, 2), 
the latter being more complex and therefore more liable to 
errors. The same dataset was then reconstructed with all 
three parameter sets. 

The results are summarized in Fig. [S] The first two 
modes are more or less restored in all cases. The third mode, 
being more complex, suffers larger distortions, and is not 
recovered with the lower inclination, where the slope of the 
phase profile is ^0.73, while the amplitudes are way off. A 
formal mode identification would give m = 1 and ^ = 3, 
totally erroneous because mode 2 already has this set of 
(^, m). There is no such problem with the higher inclination 
case, though. We conclude that although small errors in the 
binary parameters are not a serious obstacle against a suc- 
cessful EM, they may cause a misidentification of the more 
complex modes. 

Another important factor is the limb darkening. The os- 
cillation amplitude caused by a given mode in the whole-disc 
integrated flux depends significantly on the adopted limb 
darkening law, and so does the magnitude of the modula- 
tions during the eclipses. All our previous runs have been 
made with the assumption of a linear limb darkening law. 
With the full polynomial law, all modes would produce 
larger distortions during the eclipses - sometimes by orders 
of magnitude, as is the case for = 3 modes, for instance. 
The overall pattern of the distortions, however, remains the 
same. The net effect is that EM would give false amplitudes 
for the maps, but the pattern would not be modified, there- 
fore the mode identification is not hampered by small errors 
in the limb darkening. 

We did check the above conjecture by additional tests, 
although we limited ourselves to linear and polynomial limb 
darkening laws for the same stellar model, because it is 
enough to induce changes near the disc limb - the place 
where the differences really matter for the Eclipse Map- 



ping. The coefficients were taken from Claret (2001). Light 
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Figure 7. The dependence of the reconstructed profiles on the number of involved eclipses. Amplitude and phase diagrams are shown 
for each of the 3 modes. Reconstruction was made in the oblique rotator of System 1. Rows contain the following: (a) - non-resonant 
frequencies, 2 eclipse cycles; (b) - non-resonant frequencies, 10 cycles; (c) - tidally resonant frequencies, 10 cycles; (d) - uniformly splitted 
frequencies, 10 cycles. 
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Figure 8. The effect of biased geometric parameters on the reconstruction, shown in terms of amplitude and phase profiles. Each row 
shows results for a particular inclination (shown at left). Other notations are the same as in Fig. [t] 



curves were generated for a couple of selected modes = 3 
cases inclusive) with polynomial limb darkening, then eclipse 
mapped with the linear law. The other possible combina- 
tions were also performed, and the results confirm the ex- 
pected behavior. The amplitudes indeed varied greatly, but 
the nodal points of the amplitude profile and the slope of 
the phase profile both remained the same. We do not show 
them here because a similar case occurs implicitly in the 
tests of the next section. 



3.6 Hidden modes 

Due to symmetries in the pulsation patterns, for each mode 
there are certain axial orientations for which the condition 
of partial or complete cancellation occurs, that is, the in- 
tegration over the visible stellar hemisphere gives zero net 



flux variation. For example, all antisymmetric modes (with 
l — \m\ = odd) seen edge-on, as well as all sectoral modes seen 
pole-on, are subject to cancellation. For each mode there are 
also certain intermediate angles at which complete cancella- 
I Chadid et al. 12001 for an extensive study) . 



tion occurs (see 



Now in eclipsing binaries the inclination of the orbit is 
close to 90°. Antisymmetric modes on an aligned rotator 
thus will be close to the condition of complete cancellation. 
They show up during the eclipses however, because the sym- 
metry of the surface integration that led to the cancellation 
effect is lifted on a partially non-centrally occulted disc (pro- 
vided that the system is only close to, but not exactly at the 
edge-on configuration). Because it is customary in time se- 
ries analysis to ignore the data segments affected by the 
eclipses so as to avoid the complications (sidelobes and false 
peaks) caused by modulated sinusoids, such 'hidden modes' 
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Figure 9. Simulated light curves for linear (dashed lines) and 
non-linear (solid lines) limb darkening laws of the same stellar 
model. The top curves show the overall fluxes, while the oth- 
ers detail the contributions of the individual modes. Each group 
of curves has its base level marked by the dot-dashed segments 
around them. Note in particular the large differences for mode 
(3,1). 



may go undetected. Therefore their signal, amplified during 
the eclipses and unaccounted for by any frequency, certainly 
will affect the reconstruction. 

The aligned rotator with multiple modes presented in 
Section 3.3 contains in fact such a hidden mode, (3, 1) 



As 

mentioned in the previous section, = 3 modes were found 
to be extremely sensitive on the assumed limb darkening. So 
we repeated exactly that case, this time with a non-linear 
limb darkening. Fig. [9] compares the generated light curves, 
also separately showing the contributions of each individual 
mode. The second mode, (3, 1), is practically hidden, giv- 
ing almost zero net flux variation outside, but a significant 
contribution during the eclipses. 

To fully simulate the effect of hidden modes, we 
searched for frequencies in the artificial data rather than 
using the model values. For this purpose a larger dataset 
was generated, covering 5 consecutive eclipses but extend- 
ing over almost the full orbital cycle, in order to contain 
enough data for a time series analysis. We used the program 
Period04 (Lenz & Breger 2005) to derive the frequencies, 
after having cut out the segments affected by the eclipses. 
The dataset could be perfectly fitted with two frequencies 
cji = 59.1545436 and UJ2 — 65.7879857, which are very close 
to the input values (lines 1 and 3 in Table |3|. 

We then fed these two frequencies and the original 
dataset to the EM. The algorithm had trouble in achieving 
a good fit to the data. The lowest attainable values for the 
R- statistics were R ^ 128, implying = 5.6, and gave very 
messy images. Thus we had to significantly relax the fitting 
criteria from R — 1 to R — 1^0 for an acce ptab le solution; 
the also increased accordingly to 6.5. Fig. 10 shows large 
residuals that the algorithm could not fit, those being due to 
the hidden mode. That they were not falsely mapped on the 
other modes is reassuring, in fact. Indeed, the modes could 
be successfully restored, though with a higher uncertainty, 
as shown in Fig. |11| 

When we included the eclipses in the time series anal- 




Figure 11. Reconstruction of the multimode case with hidden 
mode. The top row shows the reconstructed map pairs for the 
two modes, the second and third rows contain the amplitude and 
phase proflles of the solution, and in the bottom row we show the 
maps constructed from the fltted proflles. 



ysis, even 20 frequencies were not enough to properly ac- 
count for the modulations during the eclipses. Two of these 
frequencies corresponded to modes 1 and 3, and a third fre- 
quency was near to that of the hidden mode, but its ampli- 
tude was of the same order as of other 17 peaks in the spec- 
trum. This is not unexpected, since a purely Fourier-based 
spectral analysis is only appropriate for purely harmonic 
oscillations. On the other hand, subtracting the model fit- 
ted by EM from the original dataset leaves us with a much 
cleaner residual light curve, because now EM properly ac- 
counts for the distortion effect of the eclipses on the two 
modes with known frequencies. A Fourier spectral analysis 
of these residuals is still cumbersome, with Period04 giving 
two main frequencies ^ 50.55 and 60.55 but still nothing at 
the location of the true frequency 61.547. More specific ap- 
proaches would be needed for this purpose. We mention one 
promising technique described in [Bretthorst ( 1988|, which, 
with an appropriately designed model, is expected to han- 
dle the modulations of the amplitude and the instantaneous 
frequency caused by the eclipses. 

Nevertheless, our investigation demonstrates a certain 
immunity of EM against the pollution by the hidden modes, 
in that they do not jeopardize the reconstruction of the other 
modes; moreover, EM is able to properly isolate their con- 
tribution from the detected modes. Of course, we are aware 
that the simple approach outlined above may not work as 
efficiently in real circumstances. For example, multiple hid- 
den modes may produce residuals that could be isolated 
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Figure 10. Light curves and residuals for the multimode case with hidden mode. The five echpse cycles are folded and shifted vertically 
by 0.04 from each other. The left panel shows the synthetic light curve and the fit achieved by the EM. The residuals are shown in the 
middle panel. Compare this especially with the modulations of the third mode, shown in Fig.|9] The right panel shows the residuals after 
the subtraction of the two-frequency sinusoidal fitted with Period04, for a comparison. 



only with the inclusion of an unrealistically large number 
of eclipses, or even worse, could lead to a combination that 
cannot be isolated from the signal of the detected modes. 

4 CONCLUSIONS 

In this paper we have introduced the Dynamic Eclipse Map- 
ping method, designed to reconstruct the pulsation patterns 
of non-radially oscillating components in eclipsing binaries, 
with the goal of mode identification. The method uses the 
effective surface sampling of the eclipses and provides image- 
like information on the pulsation patterns, which eventu- 
ally enables a direct mode identification, without the need 
to invoke detailed models of stellar structure and atmo- 
sphere. Only a geometric model for the binary and a few 
simple assumptions on the stellar atmosphere are needed. 
The method takes the detected frequencies and the eclipse 
light curve as input data, and furnishes pairs of images for 
each mode, completely describing its spatial and temporal 
behavior. A particular advantage of the method is that it 
can in principle operate on any wavelength range, including 
wide-band photometric data, making the datasets obtained 
in space observatories potentially useful for more than time 
series analysis. 

We have performed extensive testing of the Dynamic 
Eclipse Mapping. Based on these tests, we can make the 
following conclusions. 

(i) The reconstructions do not depend dramatically on 
the eclipse geometry. In particular, partial eclipses are not an 
obstacle, provided that their contribution is properly treated 
in terms of regularization and fitting. We believe that the use 
of an axially symmetric reference map updating scheme is 
essential in providing this rather nice property. Inclinations 
very close to 90°, however, are less favorable cases, due to 
the increasing symmetry of the eclipses' surface sampling. 

(ii) Chances for a successful mode identification decrease 
with an increasing complexity of the pulsation pattern in the 
direction perpendicular to the secondary's projected orbit, 
measured roughly by ^ — |m|. There is no such limit in the 
horizontal direction. 

(iii) Multimode pulsations can also be reconstructed, pro- 
vided that the data cover a sufficient number of individual 



eclipses. Although we only presented a case with 3 simulta- 
neous modes, there is no procedural obstacle in reconstruct- 
ing a larger number of simultaneous pulsations, if the eclipse 
coverage is large enough to allow a proper separation of all 
the detected modes. The triple mode case required 5 — 10 
eclipses, a dozen of modes would certainly need much more, 
which calls for uninterrupted space-based observations, and 
probably needs a lot more computational time as well. 

(iv) Modes resonant with the orbital motion are prob- 
lematic in that the inclusion of subsequent eclipses (after a 
certain limit, determined by the resonant cycle) will not im- 
prove the reconstruction. For single modes, one cycle may 
be enough anyway, but for multiple modes it is certainly 
a limitation that no other inversion method would be able 
to overcome. Tidally excited pulsations fall into this cate- 
gory, as well as the rotational splitting of modes in tidally 
locked systems, where the difference of the frequencies is in 
1:1 resonance with the orbital motion. 

(v) Pulsations on an oblique rotator can be reconstructed 
with similar conditions as in the aligned case, assuming that 
the direction of the rotation axis is known. Without this 
essential information, EM cannot give any useful results, 
while wrong assumptions about it will almost surely cause 
false mode identification. 

(vi) Moderate errors in the geometric parameters are tol- 
erable, although they may cause false identification of the 
more complex modes. A correct account for the limb darken- 
ing is also important. In particular, it is crucial to go beyond 
the generic linear limb darkening relation, otherwise the re- 
constructed amplitudes may become too distorted. This re- 
quirement should no cause any difficulty though, with the 
observational accuracy achievable these days. 

(vii) The method tolerates the presence of hidden modes 
and is able to properly isolate their polluting effect, so 
that they don't hinder the proper identification of the other 
modes. 

Based on the above findings, an ideal target for the method 
would be a moderately wide system, with no significant tidal 
distortions of the pulsating component (s), no resonant pul- 
sations, yet not as wide as to have an eccentric orbit, so 
that the rotation axis may be assumed aligned with the or- 
bital axis. This latter problem may be overcome by high- 
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resolution spectroscopic observations that could reveal the 
Rossiter-McLaughlin effect and infer a spin axis from it, 
as it was done for DI Herculis ( [Albrecht et al.||20Q9t. The 
effect has already been detected in RZ Cas ( Lehmann 
Mkrtichian|2008| ) , and the possibility is open for other cases 
as well. 

Our method does not require that the pulsation pat- 
terns be of spherical harmonics type; only rotational sym- 
metry must hold for the modes. In principle, modes dis- 
torted by rapid rotation of the pulsating component can 
also be investigated. Lignieres et al. ( 2006 ) have shown that 



fast rotation causes an equatorial concentration of the pulsa- 
tion amplitude, while the azimuthal structure is unchanged, 
so the modes obey the same symmetry as assumed by our 
method. In addition, it was also shown by the same authors 
that the number of horizontal surface node lines remains the 
same, with small shifts in their positions. Therefore the over- 
all topological structure of the modes is unchanged, making 
them suitable for mapping with EM. Moreover, some of the 
rotationally distorted modes have a larger disc averaging fac- 
tor than their zero rotation equivalents, making them more 
easily detectable, and ultimately allowing the detection of 
higher degree modes with ^ up to 6 — 8. Only at extremely 
high speeds does the equatorial concentration flatten the 
amplitude profile (except for a small equatorial region) to 
the extent that any topological information becomes too 
weak to recognize. We did not deal with distorted modes 
in this study, though, because we believe that the impor- 
tance of such an extension will be settled by the outcome of 
the first real- world applications. 

The simple binary model used here was appropriate 
for assessing the usability of the Dynamic Eclipse Map- 
ping method. Successful applications obviously will require 
a more detailed model for the binary, but we expect that 
the inclusion of the neglected features and effects is a com- 
putational issue, and should not endanger the success of the 
mode identification. 
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